Splitting of Majorana modes due to intervortex tunneling in a p x + ip y superconductor 
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We consider a two-dimensional (p x + ip a )-superconductor in the presence of multiple vortices, 
which support zero-energy Majorana fermion states in their cores. Intervortex tunnelings of the 
Majorana fermions lift the topological state degeneracy. Using the Bogoliubov-de Gennes equation, 
we calculate splitting of the zero-energy modes due to these tunneling events. We also discuss 
superconducting fluctuations and, in particular, their effect on the energy splitting. 
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Exotic excitations obeying non-Abelian statistics play 
a key role in topological quantum computation (TQC) 
and can emerge in a variety of condensed matter sys- 
tems 0, i i, including chiral p-wave superfluids and 
superconductors [1, Q . Such non-Abelian quasiparticles 
are believed to be realized in the A-phase of superfluid 
3 He the oxide superconductor Sr 2 Ru04 Q, p-wave 
superfluids in cold atom settings 0] as well as at the 
interfaces of s-wave superconductor and topological in- 
sulator ||. Certain vortex excitations in these systems 
support zero-energy Majorana fermions residing inside 
their cores, which lead to topological degeneracy of the 
ground states and non-Abelian statistics when multiple 
spatially separated vortices are present [H, 0] ■ 

In a spin-triplet (p x + ipj,)-superconductor, the zero- 
energy states appear in half-quantum vortices, where the 
complex phase of the condensate wavefunction winds in 
one of the spin-sectors, i.e. $(r, <p) = A(r)(p x +ip y )[\ H)+ 
e llp \ ft)]. This effectively is equivalent to a full-quantum 
vortex in a spinless (p x -Mp y )-superfluid/superconductor, 
which too may host a zero-energy state [7|. The exis- 
tence of such states has been well established through 
several approaches, including quasi-classical analysis [J], 
explicit solutions of Bogoliubov-de Gennes (BdG) equa- 



tion [0, llO, lll| and an index theorem [12j. These zero 
energy states can be occupied by Majorana fermions and 
lead to a degeneracy of the many-body ground state, 
which hinges on the exact degeneracy of the single Ma- 
jorana modes in different vortex cores. However, in the 
presence of multiple vortices, tunneling of the Majorana 
modes becomes possible and these tunneling events are 
expected to lift the ground-state degeneracy to some de- 
gree m. 

For the purposes of topological quantum computa- 
tion, it is crucial to understand the stability of Majorana 
modes against different perturbations such as fluctuation 
effects and the intervortex tunneling processes. In this 
Letter we concentrate on the latter tunneling effects in a 
two-dimensional spinless [p x + ipj,)-superconductor and 
calculate the energy splitting of the Majorana modes, in- 
cluding its dependence on the distance between vortex 
cores. We find that in addition to an exponential sup- 



pression exp(— where R and £ are the intervortex 
distance and the superconducting coherence length, re- 
spectively, the amplitude of the tunneling rate oscillates 
on the scale of the Fermi wavelength, see Eq. ([7]). These 
oscillations have important consequences for TQC which 
we discuss below. We also study thermal motion of a 
vortex which leads to smearing out the fast-oscillating 
term. We calculate the typical value of the degeneracy 
splitting (defined by its root-mean-square value) which 
determines the "decoherence" time for TQC. 

Theoretical model. Our starting point is the mean-field 
BCS Hamiltonian of spinless p x + ip y superconductor (l3j 
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which enforces the p x + ip y symmetry for the supercon- 
ducting gap: A(p) = j^(j>x + ipy)- Here z — x + iy is 
complex coordinate in the 2D plane. The BdG equation 
follows from diagonalizing Hamiltonian |T]) via Bogoli- 
ubov transformation %j}(r) = Y, n [% u n(r) + 7n u n( r )]) 
and has the form (henceforth we use H = ks = 1) 
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The Hamiltonian WedG is invariant under transforma- 
tion, cTiH-B&QOi = — 7Yg dG , which relates solutions with 
positive and negative energies [Til ] . This symmetry im- 
plies that if \E' = (u n ,v n ) T is a solution of Eq. ^ with 
eigenvalue E n , then a^* — (w*,u*) T must be a solu- 
tion with the eigenvalue (—E n ). Thus, for a zero-energy 
state, we have the constraint u = v* . In the pres- 
ence of a vortex this ensures the existence of a stable, 
symmetry-protected zero-energy state. Similar to s-wave 
superconductors 14j, the hc/2e vortex can be modeled 



2 



as A(r) = f(r)e lv , where ip is the polar angle and f(r) is 
the superconducting order-parameter profile of a vortex, 

f(r) = Ao tanh (^^j with Ao being the mean-field value 
of the superconducting order parameter. By directly 
solving BdG equation @ , one finds the energy spectrum 
for the bound states in the vortex core as E n = — uj Q n, 
where n is an integer and luq ~ Aq/ef with ep being the 
Fermi energy 16 [.The eigenstate corresponding to n = 



is given by [111 
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where vf = pF/m is the Fermi velocity, and Ji(r) is 
the Bessel function. The constant phases of u and v are 
chosen to satisfy the requirement that u* = v. Eq. (g]) 
was obtained assuming Ao <C £f which is typical for 
weak-coupling superconductors. Using the zero-mode so- 
lution (u,v) T with u = v* , we can construct the Majo- 
rana quasiparticle operator j = 7* = / d 2 r [ip(r)u*(r) + 
4>Hr)v*(r)]. 

Let us now consider the situation with 2N vortices 
pinned at positions Ri. If we ignore the fluctuation ef- 
fects and the tunneling events, the superconducting or- 
der parameter can be represented as A(r) = ]~Ii=i f( r ~ 
Ri)exp [iJ2i'Pi( r )]> where <pi(r) = arg(r - Ri). Near 
the fc-th vortex core, the phase of the order parameter can 
be approximated by (p k {r) + Q k with Q k = J2i^k <Pi(Rk)- 
Thus, one can generalize the zero-energy solution ob- 
tained for a single vortex to the situation at hand [l3[ : 
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The 2N Majorana fermions residing in the vortex cores 
can be combined in pairs to create N Dirac fermions, 
c = (ji + = (ji - ijj). This allows one 

to enumerate all degenerate ground states, which can be 
occupied by N fermions |3j. 

To understand the nature of the states created by the 
pairs of Majorana fermions, we focus here on the case 
of two vortices residing at the positions i?i and R 2 . 
Formally, a vortex potential in Eq. (J2J) is similar to a 
quantum well, in which the core excitations including 
the Majorana modes may reside. The two-vortex case 
is, therefore, similar to a double-well problem in which 
the tunneling events are expected to lead to a split- 
ting of the originally degenerate energy levels [lj|. As 
a result, the degeneracy between Dirac fermion states 
&c\l) = |1) and &c\0) = is lifted by the tunnel- 
ing between the two vortices. Here the Dirac fermion 
c = (71 + 172) /y/2 is constructed out of two Majorana 
fermions 71 and 72 . The energy difference between these 
two states is the main quantity of interest in this work. 



Similarly to a double-well problem, we first find the 
single- vortex solutions at R\ and -R2: ^1 = (ui, vi) T and 
^2 = (""2, v 2 ) T , and then construct two- vortex wavefunc- 
tions W± = (^1 ± e ta ^2) /V%, which correspond to the 
energies E+ and respectively. The energies E + and 
E_ are related by E + = —E-. The particle-hole sym- 
metry of the BdG equations requires that eri^^ = 
Combining this constraint with the properties of the zero 
energy solutions u = u*, we find that a = ir/2. Con- 
sequently, the Dirac fermion operators c and c* can be 
identified as annihilation and creation operators of the 
single-particle state i.e 
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Thus, E+ — E- = 2E+ is simply the energy splitting 
between the occupied state |1) and unoccupied state 
|0) (lo| . The energy of the state for example, can 
be calculated as the appropriate overlap integral between 
two zero-energy states [iff. When the distance between 
two vortices is much larger than superconducting coher- 
ence length R = \R 2 — R\\ ^> £ = vp/A , the wave- 
function ^2 is exponentially small in region close to the 
vortex at Ri . Then, the splitting energy E + is given by 
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Here E is the half plane x € (0, 00), y £ (—00, 00) contain- 
ing one of the vortex at R\ = (R/2,0). The other vortex 
is located at R 2 = (-R/2, 0). The integral in Eq. ([6]) can 
be calculated using the explicit form of the solution for 
^1 and We first transform integral in Eq. © over 
half-plane into a line integral over the boundary of E at 
x = 0: 
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where a = R/2£ and A = pf£,- Upon evaluating the 
integral above, we find the splitting energy to be 
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Here we neglect the corrections to the prefactor of order 
(Pf£) _1 <S 1- As one can see from Eq. 0, in addition 
to the expected exponential decay, the splitting energy 
oscillates rapidly on the Fermi wavelength. These oscil- 
lations originate from the quantum interference between 
two zero-energy eigenstates located at Ri and R 2 , see 
Fig [T] The analytical result for the degeneracy split- 
ting (|7|) is in qualitative agreement with recent numerical 
studies of the topological degeneracy in the quantum Hall 
state at Landau level filling fraction v = 5/2 171 ] as well 
as in Kitaev's honeycomb lattice model [18( ■ We note in 
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FIG. 1: (color online). Schematic plot of a two- vortex con- 
figuration. The shaded region corresponds the order param- 
eter profile f(r) with vortices located at Ri and .R2. The 
solid (blue) line represents the real part of the wavefunction 
ui(r) + iu2(r). 



passing here that the honeycomb lattice model can be 
mapped to a p-wave superconductor. However, the el- 
ementary topological excitations in the Kitacv's model 
(vortices living on a plaquette) are different from half- 
quantum vortices in the chiral p-wave superconductors 
although both have non-Abelian Ising anyons. 

In the discussion above we assumed zero temperature 
limit. At finite temperature, the error rate also comes 
from thermal fluctuations [1]. To take advantage of 
the topological quantum computation, the temperature 
should be smaller than the excitation gap, which in p- 
wave superconductors is given by the level spacing of 
bound states in the core, wo ~ A^/Ep. At T <C ljq, 
thermal population of the excited bound states in the 
vortex core is exponentially small oc exp (— ujq/T), and, 
thus, the error rate is suppressed. In addition to this well 
known argument finite temperature leads to fluctua- 
tion effects discussed in the remainder of the paper. 

Fluctuation effects. The BdG equation ^ follows from 
the mean-field Hamiltonian |T]) and assumes that the su- 
perconducting order parameter is a fixed field with no in- 
ternal dynamics. This concerns both the uniform order- 
parameter background Ao and the vortex field, which in 
Eq. @ is assumed to be an externally imposed defect 
pinned in a certain location. However, the mean-field 
Hamiltonian is an approximation and there are correc- 
tions to it due to classical and quantum fluctuations of 
the order parameter which can be separated into ampli- 
tude and phase fluctuations. The amplitude fluctuations 
manifest themselves most strongly near the transition 
and, even though they do exist in the superconducting 
phase as well, they are well-gapped and their effects are 
not as dramatic there. The bulk phase fluctuations are 
gapless in two dimensions both in a neutral system and 
in a charged supcrfluid, and, thus, can propagate over 
large distances. However, if the film thickness is finite, 
the phase fluctuations do not propagate beyond a certain 
magnetic screening length and hence become local [l9j . 



If a vortex is present, the local phase fluctuations effec- 
tively restore the vortex dynamics and, in particular, give 
rise to spatial motion of the vortex position as opposed 
to a fixed static configuration assumed in Eq. @ . This 
motion as well as any other local perturbations are un- 
likely to destroy a single topological Majorana mode, but 
the vortex motion certainly has a strong effect on the 
Majorana mode splitting, as can be seen from Eq. (J7]). 
Indeed, the splitting contains the fast-oscillating func- 
tion cos (ppR+ ?) and as such is extremely sensitive to 
the vortex positions. The Fermi wavelength is by far the 
smallest lengthscale in the problem and motion of the 
vortex effectively "smears out" the fast-oscillating term, 
introducing a new lengthscale in the problem associated 
with vortex dynamics. 

In the superfluid state at zero temperature, vortex dy- 
namics is quantum, while at higher temperatures the fluc- 
tuations of the vortex position becomes effectively classi- 
cal and have a timescale associated with it. To illustrate 
the main effect of fluctuations, we adopt here a simple 
phcnomcnological model of a single vortex moving a pin- 
ning potential, V(Ry) [13, Hl| 



1 



2m, 



{p v -A) 2 + V{R v ). 
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Here p v and Rv are the canonical momentum and co- 
ordinate of the vortex core, respectively, m v is the vor- 
tex effective mass, determined by the virtual transitions 
among the quasiparticle states caused by the vortex mo- 
tion which is treated as a phenomenological parameter 
below. Eq. JHJ also includes effective gauge field A. due 
to the Magnus force, which modifies the vortex motion, 
however we ignore it in actual calculations of fluctua- 
tion vortex dynamics. Assuming the displacement of the 
vortex position from equilibrium is small, one can ap- 
proximate the potential V(R V ) by a harmonic "pinning 
trap": V(R V ) ~ hm^LJ^Ry — R v o) 2 , where R v0 is the 
equilibrium position. The total Hamiltonian, which ac- 
counts for the coupling between vortex motion and the 
quasiparticles residing in the vortex core reads 21 1 



Tt = H. v 



Hbcs({-Rv}), 



(9) 



where Hbcs is the BCS Hamiltonian |T]). The order pa- 
rameter A(r) in Eq. (JTJ) now depends on the vortex posi- 
tions {R v }. Below, we consider the temperature regime 
lo v <C T <C u>o, where the fluctuations of the vortex po- 
sition are dominated by classical fluctuations. In this 
parameter regime the motion of the vortex is slow com- 
pared to the time scale of quasiparticle dynamics inside 
the core, lUq 1 . Therefore, the solution for the zero-energy 
eigenstates remains valid, and the calculation of the split- 
ting energy ([7|) goes through as before. However, the en- 
ergy splitting, E+ , itself becomes a random variable since 
it depends on the intervortex separation, which fluctu- 
ates. In general, we are interested in the full probability 
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distribution function of the eigenvalue splitting given by 
P[.E + ] = (5 [E+ — £? + (Hi, H 2 )]), where the average of a 
function /(Hi, H2) is defined as the integral over the vor- 
tex positions (/) = / d 2 Hid 2 H 2 /(H 1 ,H 2 )pi(Hi) /92 (H 2 ) 
weighted with the diagonal element of the density ma- 
trix of a harmonic oscillator: p(R) = exp(— H 2 /Z 2 ) /irl 2 , 
where I = tanh(/3w v /2) represents typical de- 

viation of the vortex position from the equilibrium. Be- 
low we calculate explicitly the average Majorana mode 

splitting, (£+), and its root-mean-square value, (E^_). 
In the limit lo v <C T, I is asymptotically given by 
I ~ \/2T/rn v ui^. Using the value for the vortex mass 



of Refs. [y, 21 1, we estimate that in the physically rele- 
vant regime pj, 1 <C I -C £ <C R the expression for (E^ 
becomes 
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Compared with Eq.Q, the averaged splitting rate be- 
comes exponentially smaller than its "static" value 
due to the vortex position fluctuations. This smallncss 
originates from the cosine function in Eq. ([7]) , which oscil- 
lates in sign and whose average value is therefore small. 
One may, therefore, conclude naively that the relevant 
suppression parameter is large since pfI 3> 1. This ef- 
fect is similar to the "nominal" suppression of the Friedel 
oscillations and RKKY interactions in disordered met- 
als [22[ and superconductors [23| ■ In the case of RKKY 
interaction it is well-known that the latter ensemble- 
averaged result is unphysical and the typical value of 
Friedel interaction terms is not small [22j. In our case 
the situation is similar and this is illustrated by the root- 
mean-square value 



(E* 



\/2Ao exp (— H/£) 
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which much exceeds the average (fTUj) . {E+) / J {E_ 



exp (— ^ppl 2 y This result for y {E\) represents a typical 
value of the splitting in the parameter regime of R 3> £ 
which sets the upper limit for the duration of the braiding 
operations in TQC. 

The dependence of the degeneracy splitting on physical 
parameters is important for TQC since it sheds light on 
both the adiabaticity of braiding operations and decoher- 
ence. In particular, the magnitude of the splitting energy 
sets the upper limit for the time of braiding operations 
which have to be performed adiabatically with respect 
to the excitation gap in the system u)q. On the other 
hand, because of the lifting of the topological state de- 
generacy the statistical phase gets smeared out at times 

^ <^ l/\/ Thus, the non-Abelian statistics of the 
Majorana excitations can be resolved for a wide range 



of time scales oj 1 <C t -C 1/y as l° n g as the vor- 
tices well-separted so that the splitting energy remains 
exponentially small. Furthermore, the sign of the split- 
ting energy can be either positive or negative, energeti- 
cally favoring unoccupied or occupied state by the Dirac 
fermion. Because of the oscillations of the splitting en- 
ergy on the atomic length scale, the initialization of the 
qubit in the desired quantum state as well as the read-out 
based on bringing two vortices together would become 
difficult. 

We conclude by pointing out that our work has im- 
portant consequences for topological quantum computa- 
tion using non-Abelian anyonic degeneracy as proposed 
in several related architectures recently [H, [2, H EH- 
The energy splitting we calculate will suppress topolog- 
ical immunity since it destroys the exact degeneracy of 
the quasiparticle subspace. This may adversely affect 
the fault tolerant properties of the topological quantum 
computation. Our work is the first analytical theory of 
quantum decoherence in topological quantum computa- 
tion demonstrating the limitation of the topological pro- 
tection due to quasiparticle tunneling. 
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